{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Elasticities and other comparisons\n",
    "\n",
    "In section 4.3, calculations for income and price elasticities are reported in the text. We also report the estimate of a compensating variation for being in good health. This notebook does these calculations. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib import pyplot as plt\n",
    "import numpy as np\n",
    "from scipy import stats \n",
    "from importlib import reload\n",
    "import pandas as pd\n",
    "from itertools import product\n",
    "import os \n",
    "import sys\n",
    "import pickle as pickle\n",
    "module_path = os.path.abspath(os.path.join('..'))\n",
    "if module_path not in sys.path:\n",
    "    sys.path.append(module_path)\n",
    "from model import micro, macro, params\n",
    "from model import distributions as dist"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/Users/francoislangot/opt/anaconda3/bin/python\n",
      "3.9.12 (main, Apr  5 2022, 01:53:17) \n",
      "[Clang 12.0.0 ]\n",
      "sys.version_info(major=3, minor=9, micro=12, releaselevel='final', serial=0)\n"
     ]
    }
   ],
   "source": [
    "print(sys.executable)\n",
    "print(sys.version)\n",
    "print(sys.version_info)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Loading"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "name\n",
      "sigma       2.097183\n",
      "beta        0.970000\n",
      "phi         0.304140\n",
      "psi         0.168902\n",
      "delta_h1   -0.967380\n",
      "delta_h2    3.487244\n",
      "eta         0.000000\n",
      "tfp         1.000000\n",
      "price       1.000000\n",
      "risk        0.000000\n",
      "Name: value, dtype: float64\n"
     ]
    }
   ],
   "source": [
    "df = pd.read_pickle('../estimation/output/params_ref_us.pkl')\n",
    "pars = df.loc[:,'value']\n",
    "print(pars)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Benchmark simulation for GE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sim(gr_tfp,gr_price,p,co,ge=True,rent=None,wage=None,taxrate=None):\n",
    "    outcomes = ['m','y','c','k','n','s','csy','ksy','h','g2','g3','g4',\n",
    "            'tgood','tbad','r','w','tax','oop']\n",
    "    sims = pd.Series(index=outcomes,dtype='float64')\n",
    "    ne = 10\n",
    "    nk = 100\n",
    "    nkd = 30\n",
    "    size = 2* ne * nk\n",
    "    opt = pd.DataFrame(index=np.arange(0,size),columns=['e','h','k','ps',\n",
    "                                                        'c','m','kp','v'])\n",
    "    theta = params.flexpars(sigma=p['sigma'],beta=p['beta'],\n",
    "                            phi=p['phi'],psi=p['psi'],\n",
    "                            delta_h1=p['delta_h1'],delta_h2=p['delta_h2'],eta=0.0,\n",
    "                            tfp=gr_tfp*p['tfp'],price=gr_price*p['price'])\n",
    "    # option for the numerical solution\n",
    "    m  = 2.5\n",
    "    op = params.settings(ne=ne,nk=nkd,maxk=190.0,curv=0.5,nprocs=48)\n",
    "    inc = params.incprocess(country=co)\n",
    "    inc.tauchen(ne=ne,m=m)\n",
    "    aux = params.auxpars(country=co)  \n",
    "    #Decision rules\n",
    "    if ge:\n",
    "        csumers = micro.bellman(options=op,flex=theta,aux=aux,inc=inc,rent=5.6e-2)\n",
    "    else :\n",
    "        csumers = micro.bellman(options=op,flex=theta,aux=aux,inc=inc,rent=rent,taxrate=taxrate,wage=wage)\n",
    "    csumers.compute_cash()\n",
    "    csumers.itervalue()\n",
    "    # distribution\n",
    "    stats = dist.stationary(dp=csumers,nk=nk)\n",
    "    stats.blowup()\n",
    "    stats.compute()\n",
    "    # general equilibrium\n",
    "    if ge:\n",
    "        eq = macro.equilibrium(stats=stats,taxes=False,rent=True)\n",
    "    else:\n",
    "        eq = macro.equilibrium(stats=stats,inirent=rent,taxes=False,rent=False)   \n",
    "    eq.solve()\n",
    "    aggs = eq.aggregates()\n",
    "    hlth = eq.healthreport()\n",
    "    # saving aggregate outcomes\n",
    "    sims['m']     = aggs.M\n",
    "    sims['y']     = aggs.Y\n",
    "    sims['c']     = aggs.C\n",
    "    sims['k']     = aggs.K\n",
    "    sims['n']     = aggs.N\n",
    "    sims['s']     = p['price']*aggs.M/aggs.Y\n",
    "    sims['csy']   = (aggs.C+p['price']*aggs.M)/aggs.Y\n",
    "    sims['ksy']   = aggs.K/aggs.Y\n",
    "    sims['h']     = hlth.pH\n",
    "    sims['g2']    = hlth.gradient[0]\n",
    "    sims['g3']    = hlth.gradient[1]\n",
    "    sims['g4']    = hlth.gradient[2]\n",
    "    sims['tgood'] = hlth.pTransGood\n",
    "    sims['tbad']  = hlth.pTransBad\n",
    "    sims['r']     = eq.rent\n",
    "    sims['w']     = eq.wage\n",
    "    sims['tax']   = eq.tax\n",
    "    sims['oop']   = m*p['price']*aux.copay\n",
    "    # saving decision rules\n",
    "    opt.loc[:,'ps'] = eq.stats.probs\n",
    "    for i,s in enumerate(stats.states):\n",
    "        e,h,k = s\n",
    "        opt.loc[i,['e','h','k']] = [e,h,k]  \n",
    "        opt.loc[i,'c'] = eq.stats.optc[e,h,k]\n",
    "        opt.loc[i,'m'] = eq.stats.optm[e,h,k]\n",
    "        opt.loc[i,'kp'] = eq.stats.optk[e,h,k]\n",
    "        opt.loc[i,'v'] = eq.stats.value[e,h,k]\n",
    "    return opt, sims\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "m         0.603974\n",
      "y         4.267132\n",
      "c         2.306049\n",
      "k        27.841594\n",
      "n         1.327392\n",
      "s         0.141541\n",
      "csy       0.681963\n",
      "ksy       6.524662\n",
      "h         0.879615\n",
      "g2        1.141292\n",
      "g3        1.229722\n",
      "g4        1.301956\n",
      "tgood     0.969416\n",
      "tbad      0.223414\n",
      "r         0.010194\n",
      "w         1.978646\n",
      "tax       0.202095\n",
      "oop       0.341029\n",
      "dtype: float64\n"
     ]
    }
   ],
   "source": [
    "opt_ref,sims_ref = sim(1.0,1.0,pars,'us',ge=True)\n",
    "print(sims_ref)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Value of Being in Good Health"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With a steady-state approximation of the welfare given by \n",
    "$$\n",
    "\\frac{1}{1-\\beta} \\left( \\frac{c^{1-\\sigma}}{1-\\sigma} + \\Pr(h=1) \\phi \\right),\n",
    "$$ \n",
    "we deduce that consumption reduction compensated by a rise in probability of being in good health is \n",
    "$$\n",
    "\\delta_c = 1-\\left(1-\\Delta \\Pr(h=1)\\frac{\\phi}{\\left(\\frac{c^{1-\\sigma}}{1-\\sigma}\\right)}\\right)^{1/(1-\\sigma)}\n",
    "$$ \n",
    "\n",
    "\n",
    "Given the estimated parameters $\\{\\sigma,\\phi\\}$ and the steady state values of $c$, we deduce $\\delta_c \\approx 1$\\% for $\\Delta \\Pr(h=1)= 1\\%$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.754668845032358\n"
     ]
    }
   ],
   "source": [
    "DPr_h   = 0.01\n",
    "uc      = (sims_ref['c']**(1-pars['sigma'])) / (1-pars['sigma'])\n",
    "delta_c = 1 - (1- DPr_h * pars['phi'] / uc )**(1/(1-pars['sigma']))\n",
    "print(delta_c*100)  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Partial Equilibrium"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "m         0.603813\n",
      "y         4.266882\n",
      "c         2.306383\n",
      "k        27.837339\n",
      "n         1.327392\n",
      "s         0.141512\n",
      "csy       0.682043\n",
      "ksy       6.524048\n",
      "h         0.879930\n",
      "g2        1.140729\n",
      "g3        1.228733\n",
      "g4        1.300610\n",
      "tgood     0.969416\n",
      "tbad      0.224062\n",
      "r         0.010194\n",
      "w         1.978646\n",
      "tax       0.202095\n",
      "oop       0.341029\n",
      "dtype: float64\n"
     ]
    }
   ],
   "source": [
    "opt_ref_pe,sims_ref_pe = sim(1.0,1.0,pars,'us',ge=False,rent=sims_ref['r'],wage=sims_ref['w'],taxrate=sims_ref['tax'])\n",
    "print(sims_ref_pe)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Price elasticity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "m         0.555804\n",
      "y         4.290280\n",
      "c         2.304977\n",
      "k        28.236938\n",
      "n         1.327392\n",
      "s         0.129549\n",
      "csy       0.666805\n",
      "ksy       6.581607\n",
      "h         0.879974\n",
      "g2        1.160165\n",
      "g3        1.274983\n",
      "g4        1.335073\n",
      "tgood     0.969416\n",
      "tbad      0.224148\n",
      "r         0.010194\n",
      "w         1.978646\n",
      "tax       0.202095\n",
      "oop       0.341029\n",
      "dtype: float64\n"
     ]
    }
   ],
   "source": [
    "gr_price = 1.25\n",
    "opt_rnd_pe,sims_rnd_pe = sim(1.0,gr_price,pars,'us',ge=False,rent=sims_ref['r'],wage=sims_ref['w'],taxrate=sims_ref['tax'])\n",
    "print(sims_rnd_pe)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We use an arc elasticity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-0.37261101411333836\n"
     ]
    }
   ],
   "source": [
    "eta_p = ((sims_rnd_pe['m']-sims_ref_pe['m'])/(sims_rnd_pe['m']+sims_ref_pe['m']))/((gr_price-1)/(1+gr_price))\n",
    "print(eta_p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Income Elasticity\n",
    "\n",
    "## Partial Equilibrium"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "gr_tfp = 0.75\n",
    "opt_inc_pe,sims_inc_pe = sim(gr_tfp,1.0,pars,'us',ge=False,rent=sims_ref['r'],wage=gr_tfp*sims_ref['w'],taxrate=sims_ref['tax'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.4996773903039379"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "eta_y = ((sims_inc_pe['m']-sims_ref_pe['m'])/(sims_inc_pe['m']+sims_ref_pe['m']))/((gr_tfp-1)/(1+gr_tfp))\n",
    "eta_y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## General Equilibrium"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "opt_inc_ge,sims_inc_ge = sim(gr_tfp,1.0,pars,'us',ge=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.5729752584198164"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "eta_y = ((sims_inc_ge['m']-sims_ref['m'])/(sims_inc_ge['m']+sims_ref['m']))/((gr_tfp-1)/(1+gr_tfp))\n",
    "eta_y"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
